phosphoDataset msqrobPTM final workflow with GP

library(QFeatures)
library(msqrob2)
library(tibble)
library(tidyr)
library(dplyr)
library(gt)
library(plotly)
library(stageR)
library(poolr)
library(RColorBrewer)
library(seqinr)
library(stringr)
library(ExploreModelMatrix)
library(data.table)

1 Import data

df <- read.csv("evidence.txt",header=TRUE, sep = ",")
#I can ignore the sigma proteins (those were just used for QC purposes)
df <- df %>% filter(!grepl("Sigma", df$Proteins))
#filter out phospho sites with less than 75% probability
filter_rows <- sapply(df$Phospho..STY..Probabilities, function(x){
  if (!x==""){
    prob <- str_extract_all(x,  "(?<=\\().+?(?=\\))")[[1]]
    if (!any(as.double(prob) > 0.75)){return(which(df$Phospho..STY..Probabilities == x))}
  }
}, USE.NAMES = F)
filter_rows <- unique(unlist(filter_rows))
nrow(df)

## [1] 178235

df <- df[-filter_rows,]
nrow(df)

## [1] 171780

#Now get format into wide format
df_wide <- pivot_wider(df, id_cols = c("Sequence", "Modifications", "Modified.sequence", "Proteins", "Leading.proteins",
                                        "Reverse", "Potential.contaminant", "Protein.group.IDs", "Leading.razor.protein"),
                       names_from = "Raw.file", names_prefix = "Intensity_", values_from = "Intensity", values_fn = max)

ecols <- grep("Intensity",colnames(df_wide))
#order dataframe by protein, for the normalisation step
df_wide = df_wide[order(df_wide$Leading.razor.protein),]
pe <- readQFeatures(df_wide,ecol= ecols,name="peptidoformRaw")

rownames(pe[["peptidoformRaw"]]) <- rowData(pe[["peptidoformRaw"]])$Modified.sequence

1.1 Experimental Layout

metadata <- read.csv("Experimental Design.csv")

colData(pe)$file <- sapply(str_split(rownames(colData(pe)), "_"), function(x) x[[2]])
metadata$File <- sub(x=metadata$File, pattern=".raw",replacement = "")
colData(pe)$condition <- sapply(colData(pe)$file, function(x){
  metadata[metadata$File==x,]$Condition
}, USE.NAMES = F)
colData(pe)$subset <- sapply(colData(pe)$file, function(x){
  metadata[metadata$File==x,]$Subset
}, USE.NAMES = F)
colData(pe)$condition <- case_when(grepl("A", colData(pe)$condition) ~ "A",
                                        TRUE ~ "B")
colData(pe)

## DataFrame with 90 rows and 3 columns
##                            file   condition      subset
##                     <character> <character> <character>
## Intensity_QX21595GM   QX21595GM           A           y
## Intensity_QX21643GM   QX21643GM           A           y
## Intensity_QX21469GM   QX21469GM           B           x
## Intensity_QX21484GM   QX21484GM           B           x
## Intensity_QX21487GM   QX21487GM           A           x
## ...                         ...         ...         ...
## Intensity_QX21782GM   QX21782GM           A           y
## Intensity_QX21794GM   QX21794GM           A           y
## Intensity_QX21803GM   QX21803GM           B           y
## Intensity_QX21806GM   QX21806GM           A           y
## Intensity_QX21812GM   QX21812GM           A           y

1.2 Protein dataset

df_prot <- read.csv("Non-enriched/evidence.txt",header=TRUE, sep = ",")
#I can ignore the pool samples (those were just used for QC purposes)
#2 samples lack metadata info, so I will remove them for now as well
df_prot <- df_prot %>% filter(!grepl("Sigma", df_prot$Proteins),
                    #extra files not present in the ptm dataset, so not useful
                    !grepl("QE37656GM|QE37749GM|QE37769GM|QE37811GM|QE37913GM|QE37698GM|QE37781GM|QE37919GM|QE37898GM|QE37979GM|QE37713GM", df_prot$Raw.file))

#Now get format into wide format
df_prot_wide <- pivot_wider(df_prot, id_cols = c("Sequence", "Modifications", "Modified.sequence", "Proteins", "Leading.proteins",
                                        "Reverse", "Potential.contaminant", "Protein.group.IDs", "Leading.razor.protein"),
                       names_from = "Raw.file", names_prefix = "Intensity_", values_from = "Intensity", values_fn = max)

ecols <- grep("Intensity",colnames(df_prot_wide))
#order dataframe by protein, for the normalisation step
df_prot_wide = df_prot_wide[order(df_prot_wide$Leading.razor.protein),]
prot <- readQFeatures(df_prot_wide,ecol= ecols,name="proteinRaw")
rownames(prot[["proteinRaw"]]) <- rowData(prot[["proteinRaw"]])$Modified.sequence

Total protein overlap

prot_GP <- unique(df_prot_wide$Leading.razor.protein)
prot_P <- unique(df_wide$Leading.razor.protein)
#how many of the protein in enriched dataset are not in non-enriched dataset?
n_diff <- length(setdiff(prot_P, prot_GP))
n_diff/length(prot_P)

## [1] 0.2466443

25% of proteins in the enriched dataset do not find a counterpart in the non-enriched dataset

1.2.1 Experimental Layout

metadata <- read.csv("Non-enriched/metadata.csv")

colData(prot)$file <- sapply(str_split(rownames(colData(prot)), "_"), function(x) x[[2]])
metadata$File <- sapply(strsplit(metadata$File, ".", fixed = T), function(x) x[[1]])
metadata <- metadata %>% filter(!grepl("QE37656GM|QE37749GM|QE37769GM|QE37811GM|QE37913GM|QE37698GM|QE37781GM|QE37919GM|QE37898GM|QE37979GM|QE37713GM", metadata$File))
colData(prot)$condition <- sapply(colData(prot)$file, function(x){
  metadata[metadata$File==x,]$Condition
}, USE.NAMES = F)
colData(prot)$subset <- sapply(colData(prot)$file, function(x){
  metadata[metadata$File==x,]$Subset
}, USE.NAMES = F)
colData(prot)$condition <- case_when(grepl("A", colData(prot)$condition) ~ "A",
                                        TRUE ~ "B")
colData(prot)

## DataFrame with 90 rows and 3 columns
##                            file   condition      subset
##                     <character> <character> <character>
## Intensity_QE37737GM   QE37737GM           A           y
## Intensity_QE37799GM   QE37799GM           B           x
## Intensity_QE37829GM   QE37829GM           A           x
## Intensity_QE37886GM   QE37886GM           A           y
## Intensity_QE37895GM   QE37895GM           A           y
## ...                         ...         ...         ...
## Intensity_QE37988GM   QE37988GM           B           x
## Intensity_QE37892GM   QE37892GM           B           y
## Intensity_QE37907GM   QE37907GM           B           y
## Intensity_QE37991GM   QE37991GM           A           y
## Intensity_QE37949GM   QE37949GM           B           x

2 Preprocessing

2.1 PTM dataset

rowData(pe[["peptidoformRaw"]])$nNonZero <- rowSums(assay(pe[["peptidoformRaw"]]) > 0, na.rm = T)
pe <- zeroIsNA(pe, i = "peptidoformRaw")
pe <- filterFeatures(pe, ~Reverse != "+")
pe <- filterFeatures(pe, ~Potential.contaminant != "+")
pe <- logTransform(pe, base = 2, i = "peptidoformRaw", name = "peptidoformLog")
pe <- pe[rowData(pe[["peptidoformRaw"]])$nNonZero>2,,]
pe <- QFeatures::normalize(pe, method = "center.median", i = "peptidoformLog", name = "peptidoform")
colData(pe[["peptidoform"]]) <- colData(pe)

2.2 Protein dataset

rowData(prot[["proteinRaw"]])$nNonZero <- rowSums(assay(prot[["proteinRaw"]]) > 0, na.rm = T)
prot <- zeroIsNA(prot, i = "proteinRaw")
prot <- logTransform(prot, base = 2, i = "proteinRaw", name = "proteinLog")
prot <- prot[rowData(prot[["proteinRaw"]])$nNonZero>2,,]
prot <- QFeatures::normalize(prot, method = "center.median", i = "proteinLog", name = "protein")
colData(prot[["protein"]]) <- colData(prot)

2.3 Normalisation via robust summarisation

prot <- aggregateFeatures(prot,
 i = "protein",
 fcol = "Leading.razor.protein",
 na.rm = TRUE,
 name = "proteinRobust",
 fun = MsCoreUtils::robustSummary)

## Your quantitative data contain missing values. Please read the relevant
## section(s) in the aggregateFeatures manual page regarding the effects
## of missing values on data aggregation.

I have to change the colnames of the proteinRobust assay so that it can be used as a centering for the ptm dataset. I have matched the corresponding samples to each other via the metadata files (correspondance/condition column) and just changed the raw files names to the raw files names of the ptm dataset. See also rawfilenames_change.csv

rawfilenames <- read.csv("rawfilenames_change.csv")
#get current protein colnames without "Intensity_"
current_colnames <- sapply(strsplit(colnames(assay(prot[["proteinRobust"]])), "_"), function(x) x[[2]])
#replace with corresponding ptm colnames
colnames(assay(prot[["proteinRobust"]], withDimnames = F)) <- 
  sapply(current_colnames, function(x){
                  ptm_colname <- rawfilenames[rawfilenames[["protein.dataset"]]==x,]$`ptm.dataset`
                  paste0("Intensity_", ptm_colname)
                }, USE.NAMES = F)

#Only take pepforms that have a parent protein
#(when there is no global profiling dataset, this will be all peptidoforms)
pepWithProtein <- rowData(pe[["peptidoform"]])$Proteins %in% rownames(prot[["proteinRobust"]])
pePepWithProtein  <- pe[["peptidoform"]][pepWithProtein,]
pe <- addAssay(pe,pePepWithProtein,"pepformRel")
#normalisation for protein abundance step
assay(pe[["pepformRel"]]) <- assay(pe[["pepformRel"]]) - assay(prot[["proteinRobust"]], 
                                                               withDimnames = F)[rowData(pe[["pepformRel"]])$Proteins,colnames(assay(pe[["pepformRel"]]))]
boxplot(assay(pe[["pepformRel"]]))

3 Differential peptidoform usage (DPFU)

3.1 Hypothesistest for each contrast

condition A vs B

colData(pe)$condition <- relevel(as.factor(colData(pe)$condition), ref = "B")
colData(pe)$subset <- as.factor(colData(pe)$subset)
colData(pe[["pepformRel"]]) <- colData(pe)

pe <- msqrob2::msqrob(object = pe, i = "pepformRel", formula = ~condition*subset, robust=FALSE)
rowData(pe[["pepformRel"]])$msqrobModels[[2]] %>%
                  getCoef

##        (Intercept)         conditionA            subsety conditionA:subsety 
##          1.6565104         -0.5246633          0.8686482          1.2166721

fraction_of_ys <- (colData(pe) %>% as_tibble() %>% filter(subset == "y") %>% nrow()) / nrow(colData(pe))
contrasts <- c("conditionA", "conditionA + conditionA:subsety", "conditionA:subsety", "conditionA + 0.5 * conditionA:subsety", "conditionA + 0.6666667 * conditionA:subsety")
L <- makeContrast(c("conditionA = 0",
                    "conditionA + conditionA:subsety = 0",
                    "conditionA:subsety = 0",
                    "conditionA + 0.5 * conditionA:subsety = 0",
                    "conditionA + 0.6666667 * conditionA:subsety = 0"),
                  parameterNames = rowData(pe[["pepformRel"]])$msqrobModels[[2]] %>%
                  getCoef %>%
                  names)
pe <- hypothesisTest(object = pe, i = "pepformRel", contrast = L)

3.1.1 Volcanoplot

volcano <- list()
for (contrast in contrasts){
volcano[[contrast]] <- rowData(pe[["pepformRel"]])[[contrast]]%>%
            ggplot(aes(x = logFC, y = -log10(pval), 
                   color = adjPval < 0.05,
                   annotation=rowData(pe[["pepformRel"]])[,3])) +
  geom_point(cex = 2.5) +
  scale_color_manual(values = alpha(c("black", "red"), 0.5)) + 
  theme_minimal() +
  ylab("-log10(pvalue)") +
  ggtitle(contrast)
}
volcano

## $conditionA

## Warning: Removed 665 rows containing missing values (geom_point).

## 
## $`conditionA + conditionA:subsety`

## Warning: Removed 665 rows containing missing values (geom_point).

## 
## $`conditionA:subsety`

## Warning: Removed 665 rows containing missing values (geom_point).

## 
## $`conditionA + 0.5 * conditionA:subsety`

## Warning: Removed 665 rows containing missing values (geom_point).

## 
## $`conditionA + 0.6666667 * conditionA:subsety`

## Warning: Removed 665 rows containing missing values (geom_point).

3.1.2 Significant peptidoforms for each contrast

tables <- list()
for (contrast in contrasts){
sigTable <- rowData(pe[["pepformRel"]])[[contrast]]
if(nrow(sigTable <- sigTable %>%
  na.exclude %>%
  filter(adjPval<0.05))>0){
sigTable <- sigTable %>%
  na.exclude %>%
  filter(adjPval<0.05) %>%
  arrange(pval) %>%
  mutate(
    se = format(se, digits = 2), 
    df = format(df, digits =2),
    t = format(t, digits = 2),
    adjPval = format(adjPval, digits = 3),
    rank = 1:length(logFC) 
  ) 
sigTable_print <- sigTable %>% mutate(peptidoform = rownames(sigTable)) %>% gt() %>% tab_header(title = md(contrast))
tables[[contrast]] <- sigTable
#knitr::kable(sigTable, caption = contrast)
}
}
knitr::kable(tables)

logFC se df t pval adjPval rank LVGGPM(Oxidation (M))DAS(Phospho (STY))VEEEGVRR -1.2542955 0.25 87 -5.1 0.0000023 0.00249 1 VS(Phospho (STY))REFHSHEFHSHEDM(Oxidation (M))LVVDPK -2.0161711 0.40 76 -5.0 0.0000035 0.00249 2 EFHSHEFHS(Phospho (STY))HEDM(Oxidation (M))LVVDPK -1.8682595 0.41 87 -4.5 0.0000194 0.00913 3 LPIVNFDYS(Phospho (STY))M(Oxidation (M))EEK -1.0615646 0.24 80 -4.4 0.0000393 0.01390 4 GHPQEESEESNVS(Phospho (STY))MASLGEK -1.1730632 0.27 49 -4.4 0.0000618 0.01748 5 LVGGPMDAS(Phospho (STY))VEEEGVRR -1.4414466 0.35 90 -4.1 0.0000965 0.02274 6 S(Phospho (STY))GEATDGARPQALPEPMQESK -1.3518899 0.34 85 -4.0 0.0001483 0.02798 7 EIPAWVPFDPAAQITK -2.0169090 0.45 25 -4.4 0.0001583 0.02798 8 EFHS(Phospho (STY))HEFHS(Phospho (STY))HEDM(Oxidation (M))LVVDPK -1.6948425 0.44 81 -3.9 0.0002024 0.03180 9 GILAADESTGS(Phospho (STY))IAKR -0.7129553 0.19 81 -3.8 0.0002476 0.03344 10 EFHSHEFHS(Phospho (STY))HEDMLVVDPK -2.1531285 0.56 81 -3.8 0.0002602 0.03344 11

logFC se df t pval adjPval rank LVGGPM(Oxidation (M))DAS(Phospho (STY))VEEEGVRR -0.9202058 0.20 87 -4.6 0.0000136 0.0193 1 VS(Phospho (STY))REFHSHEFHSHEDM(Oxidation (M))LVVDPK -1.4202699 0.32 76 -4.4 0.0000358 0.0220 2 LPIVNFDYS(Phospho (STY))M(Oxidation (M))EEK -0.8670566 0.20 80 -4.3 0.0000466 0.0220 3 GHPQEESEESNVS(Phospho (STY))MASLGEK -0.9505352 0.22 49 -4.3 0.0000690 0.0244 4 GILAADESTGS(Phospho (STY))IAK -0.7224652 0.18 73 -4.0 0.0001662 0.0470 5

4 DPTM

4.1 ptm summarisatie

Summarise peptidoforms to ptm level

4.2 Determine location of the modification

fasta <- "human.fasta"
parsed_fasta <- read.fasta(file = fasta, seqtype = "AA", as.string = T)

#modified sequence column contains _ that does not matter, but hinders location determination
rowData(pe[["pepformRel"]])$modified_sequence <- gsub("_", "", rowData(pe[["pepformRel"]])$Modified.sequence)

get_ptm_location <- function(feature, data, fasta, mod_column = "Modifications", 
                             peptide_seq_column = "Sequence", mod_seq_column = "modified_sequence", 
                             protein_column = "Leading.razor.protein", split = ",", collapse = ", "){
  prot <- data[feature,][[protein_column]]
  pep_seq <- data[feature,][[peptide_seq_column]]
  mod_seq <- data[feature,][[mod_seq_column]]
  prot_seq <- fasta[[prot]][1]

  #find location of peptide in protein
  #add -1 here, so that the addition of the location later on is correct
  pep_location <- unlist(lapply(gregexpr(pattern = pep_seq, prot_seq), min)) - 1
  final_mod <- c()
  j <- mod_seq
  #go over the modifications inside the modified sequence
  for(mod in regmatches(mod_seq, gregexpr("\\(.*?\\)\\)", mod_seq, perl=T))[[1]]){
    #find location of modification in peptide
    mod_location <- unlist(lapply(gregexpr(pattern = mod, j, fixed = T), min))
    #get location in protein (-1 because else it gives you the location after because of the presence of the modification in the string)
    location <- mod_location + pep_location -1
    #add location to modification
    mod_ <- paste(mod, location)
    #save modification
    final_mod <- c(final_mod, mod_)
    #now remove current modification from the sequence, so that we can continue to the next mod
    str_sub(j, mod_location, nchar(mod)+mod_location-1) <- ""
  }
  return(paste(final_mod, collapse = collapse))
}

rowData(pe[["pepformRel"]])$mod <- sapply(rownames(rowData(pe[["pepformRel"]])), function(x){
    get_ptm_location(x, rowData(pe[["pepformRel"]]), parsed_fasta)
}, USE.NAMES = F)

#Add ptm variable = protein + modification
rowData(pe[["pepformRel"]])$ptm <- ifelse(rowData(pe[["pepformRel"]])$mod != "",
                                               paste(rowData(pe[["pepformRel"]])$Leading.razor.protein, rowData(pe[["pepformRel"]])$mod, sep="_"), 
                                               "")

4.2.1 Get ptm level intensity matrix

Get matrix for DPTMx -> get summarised intensity values

prots <- unique(rowData(pe[["pepformRel"]])$Leading.razor.protein)
#Do for each protein
ptms <- sapply(prots, function(i) {
  pe_sub <- pe[["pepformRel"]][grepl(i, rowData(pe[["pepformRel"]])$Leading.razor.protein, fixed = T),]
  #pe_sub <- filterFeatures(pe,~grepl(Leading.razor.protein,pattern=i,fixed = T))
  #Get all unique modification present on that protein
  mods <- unique(unlist(strsplit(rowData(pe_sub)$mod, split = ", ", fixed = TRUE)))
  #Add protein info to mods
  ptm <- paste(rep(i, length(mods)), mods)
  #return all the protein-mods combinations
  ptm
})
ptms <- as.vector(unlist(ptms))

#For each ptm do
ptm_x_assay <- sapply(seq(1:length(ptms)), function(i){ 
  x <- ptms[i]
  #Get current protein and mod from ptm
  prot <- str_split(x, " ", 2)[[1]][1]
  current_ptm <- str_split(x, " ", 2)[[1]][2]
  #filter on that protein and on that mod to obtain all peptidoforms that correspond to the ptm
  pe_sub <- pe[["pepformRel"]][grepl(prot, rowData(pe[["pepformRel"]])$Leading.razor.protein, fixed = T),]
  #pe_sub <- filterFeatures(pe,~grepl(Leading.razor.protein,pattern=prot, fixed = T))
  ptm_sub <- pe_sub[grepl(current_ptm, rowData(pe_sub)$mod, fixed = T),]
  #ptm_sub <- filterFeatures(pe_sub,~grepl(mod,pattern=current_ptm, fixed=T))[["peptidoformNorm"]]
  #Get intensity values of those peptidoforms
  y <- assay(ptm_sub)
  #And summarise them to 1 row of intensity values: 1 value per sample for that ptm
  ptm_y <- try(MsCoreUtils::robustSummary(y), silent = T)
  if (is(ptm_y, "try-error")){
    ptm_y <- rep(NA, ncol(y))}
  ptm_y
})

## Warning in rlm.default(X, expression, ...): 'rlm' failed to converge in 20 steps

## Warning in rlm.default(X, expression, ...): 'rlm' failed to converge in 20 steps

## Warning in rlm.default(X, expression, ...): 'rlm' failed to converge in 20 steps

## Warning in rlm.default(X, expression, ...): 'rlm' failed to converge in 20 steps

#Then we get the intensity assay on ptm level
ptm_x_assay <- t(ptm_x_assay)
rownames(ptm_x_assay) <- ptms

4.2.2 Add to QFeatures object

Filter out ptms with all zero intensities

print(paste(nrow(ptm_x_assay), "ptms before filtering"))

## [1] "430 ptms before filtering"

filtering <- rowSums(ptm_x_assay != 0, na.rm=TRUE) > 0 
ptm_x_assay <- ptm_x_assay[filtering,]
print(paste(nrow(ptm_x_assay), "ptms after filtering"))

## [1] "427 ptms after filtering"

all(rownames(colData(pe)) == colnames(ptm_x_assay))

## [1] TRUE

rowdata <- data.frame(ptm = rownames(ptm_x_assay))
rowdata$protein <- sapply(str_split(rowdata$ptm, pattern=" "),function(x) x[1])
ptm_y_assay <- SummarizedExperiment(assays=as.matrix(ptm_x_assay), rowData=rowdata, colData=colData(pe))

pe <- addAssay(pe, ptm_y_assay, "ptmRel")

4.3 Hypothesistest for each contrast

colData(pe[["ptmRel"]])

## DataFrame with 90 rows and 3 columns
##                            file condition   subset
##                     <character>  <factor> <factor>
## Intensity_QX21595GM   QX21595GM         A        y
## Intensity_QX21643GM   QX21643GM         A        y
## Intensity_QX21469GM   QX21469GM         B        x
## Intensity_QX21484GM   QX21484GM         B        x
## Intensity_QX21487GM   QX21487GM         A        x
## ...                         ...       ...      ...
## Intensity_QX21782GM   QX21782GM         A        y
## Intensity_QX21794GM   QX21794GM         A        y
## Intensity_QX21803GM   QX21803GM         B        y
## Intensity_QX21806GM   QX21806GM         A        y
## Intensity_QX21812GM   QX21812GM         A        y

pe <- msqrob2::msqrob(object = pe, i = "ptmRel", formula = ~condition*subset, robust=FALSE, overwrite = T)
rowData(pe[["ptmRel"]])$msqrobModels[[2]] %>%
                  getCoef

##        (Intercept)         conditionA            subsety conditionA:subsety 
##         -0.5032653         -0.1630285          1.2483607         -0.2654357

contrasts <- c("conditionA", "conditionA + conditionA:subsety", "conditionA:subsety", "conditionA + 0.5 * conditionA:subsety", "conditionA + 0.6666667 * conditionA:subsety")
L <- makeContrast(c("conditionA = 0",
                    "conditionA + conditionA:subsety = 0",
                    "conditionA:subsety = 0",
                    "conditionA + 0.5 * conditionA:subsety = 0",
                    "conditionA + 0.6666667 * conditionA:subsety = 0"),   
                  parameterNames = rowData(pe[["ptmRel"]])$msqrobModels[[3]] %>%
                  getCoef %>%
                  names)
pe <- hypothesisTest(object = pe, i = "ptmRel", contrast = L, overwrite = T)

4.3.1 Volcanoplot

volcanos <- list()
for (contrast in contrasts){
library(plotly)
volcanos[[contrast]] <- 
  rowData(pe[["ptmRel"]])[[contrast]]%>%
  ggplot(aes(x = logFC, y = -log10(pval), 
             color = adjPval < 0.05,
             annotation=rowData(pe[["ptmRel"]])[,3])) +
  geom_point(cex = 2.5) +
  scale_color_manual(values = alpha(c("black", "red"), 0.5)) + theme_minimal() +
  ylab("-log10(pvalue)") +
  ggtitle(contrast)
}
volcanos

## $conditionA

## Warning: Removed 57 rows containing missing values (geom_point).

## 
## $`conditionA + conditionA:subsety`

## Warning: Removed 57 rows containing missing values (geom_point).

## 
## $`conditionA:subsety`

## Warning: Removed 57 rows containing missing values (geom_point).

## 
## $`conditionA + 0.5 * conditionA:subsety`

## Warning: Removed 57 rows containing missing values (geom_point).

## 
## $`conditionA + 0.6666667 * conditionA:subsety`

## Warning: Removed 57 rows containing missing values (geom_point).

4.3.2 Significant ptms for each contrast

tables <- list()
for (contrast in contrasts){
sigTable <- rowData(pe[["ptmRel"]])[[contrast]]
if(nrow(sigTable %>%
  na.exclude %>%
  filter(adjPval<0.05)) > 0){
sigTable <- sigTable %>%
  na.exclude %>%
  filter(adjPval<0.05) %>%
  arrange(pval) %>%
  mutate(
    se = format(se, digits = 2), 
    df = format(df, digits =2),
    t = format(t, digits = 2),
    adjPval = format(adjPval, digits = 3),
    rank = 1:length(logFC) 
  ) 
sigTable_print <- sigTable %>% mutate(peptidoform = rownames(sigTable)) %>% gt() %>% tab_header(title = md(contrast))
tables[[contrast]] <- sigTable
}
}
knitr::kable(tables)

logFC se df t pval adjPval rank sp|P10451|OSTP_HUMAN (Oxidation (M)) 284 -1.4659656 0.31 89 -4.8 0.0000075 0.00279 1 sp|P01034|CYTC_HUMAN (Oxidation (M)) 40 -0.9714811 0.22 88 -4.3 0.0000393 0.00509 2 sp|O94769|ECM2_HUMAN (Oxidation (M)) 76 -1.0615646 0.24 79 -4.3 0.0000413 0.00509 3 sp|P10645|CMGA_HUMAN (Phospho (STY)) 161 -1.2783351 0.30 61 -4.2 0.0000795 0.00735 4 sp|P10645|CMGA_HUMAN (Phospho (STY)) 142 -1.2574547 0.31 84 -4.1 0.0001125 0.00753 5 sp|P04075|ALDOA_HUMAN (Phospho (STY)) 39 -0.7100813 0.18 88 -4.0 0.0001222 0.00753 6 sp|P05060|SCG1_HUMAN (Phospho (STY)) 317 -1.1385927 0.30 70 -3.8 0.0002791 0.01475 7 sp|P10451|OSTP_HUMAN (Phospho (STY)) 280 -1.7898662 0.50 89 -3.6 0.0005233 0.02275 8 sp|O94769|ECM2_HUMAN (Phospho (STY)) 75 -1.1008025 0.31 89 -3.6 0.0005534 0.02275 9 sp|P02774|VTDB_HUMAN (Phospho (STY)) 95 -0.7719926 0.22 63 -3.6 0.0006759 0.02501 10 sp|P01042|KNG1_HUMAN (Phospho (STY)) 275 -1.3530778 0.38 39 -3.6 0.0009266 0.02893 11 sp|P51693|APLP1_HUMAN (Phospho (STY)) 515 -1.2953037 0.38 83 -3.4 0.0009384 0.02893 12 sp|P10451|OSTP_HUMAN (Phospho (STY)) 27 -1.2549447 0.37 89 -3.4 0.0010416 0.02964 13 sp|P10645|CMGA_HUMAN (Phospho (STY)) 218 -0.9143803 0.27 89 -3.4 0.0011272 0.02978 14 sp|P13521|SCG2_HUMAN (Phospho (STY)) 106 -1.2377539 0.36 48 -3.4 0.0012072 0.02978 15 sp|P02765|FETUA_HUMAN (Oxidation (M)) 321 -1.4433264 0.45 79 -3.2 0.0017541 0.03823 16 sp|P09972|ALDOC_HUMAN (Phospho (STY)) 39 -1.1682749 0.36 73 -3.2 0.0019930 0.03823 17 sp|P19823|ITIH2_HUMAN (Oxidation (M)) 64 -2.4624647 0.69 19 -3.6 0.0020215 0.03823 18 sp|P10451|OSTP_HUMAN (Phospho (STY)) 24 -0.9663830 0.30 62 -3.2 0.0021278 0.03823 19 sp|P10909|CLUS_HUMAN (Phospho (STY)) 210 -1.3693705 0.42 38 -3.3 0.0021492 0.03823 20 sp|P00747|PLMN_HUMAN (Phospho (STY)) 358 -0.9781689 0.30 58 -3.2 0.0021699 0.03823 21 sp|P24592|IBP6_HUMAN (Phospho (STY)) 152 -0.8601732 0.28 86 -3.1 0.0025077 0.04098 22 sp|P10451|OSTP_HUMAN (Phospho (STY)) 195 -1.2466403 0.40 76 -3.1 0.0025475 0.04098 23 sp|P01034|CYTC_HUMAN (Phospho (STY)) 43 -0.9600921 0.31 90 -3.1 0.0027547 0.04247 24 sp|P24593|IBP5_HUMAN (Phospho (STY)) 124 -0.8831684 0.29 66 -3.1 0.0029382 0.04348 25 sp|P10645|CMGA_HUMAN (Phospho (STY)) 370 -0.9174091 0.29 35 -3.2 0.0032049 0.04561 26 sp|P61769|B2MG_HUMAN (Phospho (STY)) 108 -0.6020597 0.20 87 -3.0 0.0034231 0.04592 27 sp|P30086|PEBP1_HUMAN (Phospho (STY)) 52 -0.9909364 0.33 80 -3.0 0.0034749 0.04592 28 sp|P04075|ALDOA_HUMAN (Phospho (STY)) 36 -1.2092350 0.40 79 -3.0 0.0037198 0.04746 29 sp|P05060|SCG1_HUMAN (Phospho (STY)) 311 -0.8607261 0.29 90 -3.0 0.0039602 0.04884 30

logFC se df t pval adjPval rank sp|P04075|ALDOA_HUMAN (Phospho (STY)) 39 -0.6427094 0.14 88 -4.5 0.0000246 0.00905 1 sp|O94769|ECM2_HUMAN (Oxidation (M)) 76 -0.8670566 0.20 79 -4.3 0.0000489 0.00905 2 sp|P01034|CYTC_HUMAN (Oxidation (M)) 40 -0.7164067 0.18 88 -3.9 0.0001623 0.02002 3 sp|P10451|OSTP_HUMAN (Oxidation (M)) 284 -0.9637083 0.25 89 -3.8 0.0002268 0.02098 4 sp|P10645|CMGA_HUMAN (Phospho (STY)) 142 -0.9261445 0.25 84 -3.7 0.0003919 0.02900 5 sp|P02774|VTDB_HUMAN (Phospho (STY)) 95 -0.6487798 0.18 63 -3.7 0.0005014 0.03092 6 sp|P30086|PEBP1_HUMAN (Phospho (STY)) 52 -0.9408715 0.26 80 -3.6 0.0006198 0.03276 7 sp|P09972|ALDOC_HUMAN (Phospho (STY)) 39 -1.0004081 0.29 73 -3.5 0.0008802 0.04071 8 sp|Q14515|SPRL1_HUMAN (Oxidation (M)) 276 -1.5307028 0.43 35 -3.5 0.0012155 0.04822 9 sp|P10645|CMGA_HUMAN (Phospho (STY)) 218 -0.7337973 0.22 89 -3.3 0.0013032 0.04822 10 sp|P10645|CMGA_HUMAN (Phospho (STY)) 161 -0.8365480 0.25 61 -3.3 0.0014779 0.04971 11 sp|P04075|ALDOA_HUMAN (Phospho (STY)) 36 -1.0796746 0.33 79 -3.3 0.0016171 0.04986 12

5 PLOTS

5.0.1 Lineplots

5.0.1.1 PTM - level

ptm_list1 <- c()
for (contrast in contrasts){
sigTable <- rowData(pe[["ptmRel"]])[[contrast]] %>% filter(adjPval<0.05)
ptm_list1 <- c(ptm_list1, rownames(sigTable))
}
ptm_list1 <- unique(ptm_list1)

plots <- list()
for (i in ptm_list1){
prot_ = str_split(i, " ")[[1]][1]
site_ = str_split_fixed(i, " ", 2)[,2]

pepform_df <- longFormat(pe[,,"pepformRel"], rowvars = c("Leading.razor.protein", "ptm"), colvars = c("condition", "subset")) %>% as.data.frame()
pepform_df <- pepform_df %>% filter(Leading.razor.protein == prot_)
pepform_df <- pepform_df %>% filter(grepl(site_, ptm, fixed = T))
pepform_df$FeatureType <- "Peptide"
pepform_df <- pepform_df %>% arrange(condition, subset, rowname)

ptm_df <- longFormat(pe[,,"ptmRel"], rowvars = c("protein", "ptm"), colvars = c("condition", "subset")) %>% as.data.frame()
ptm_df <- ptm_df %>% filter(ptm == i)
ptm_df$FeatureType <- "PTM"
ptm_df <- ptm_df %>% arrange(condition, subset, rowname)

ptm_estimate <- ptm_df %>% select(c("primary", "colname", "condition", "subset")) %>%
                           mutate(rowname = paste(ptm_df$rowname, "estimate"),
                                  ptm = paste(ptm_df$ptm, "estimate"),
                                  Protein = paste(ptm_df$Protein, "estimate"),
                                  FeatureType = "PTM_estimated",
                                  value = NA,
                                  assay = "model")
fixeffects <- rowData(pe[["ptmRel"]])$msqrobModels[[unique(ptm_df$ptm)]] %>% getCoef
ptm_estimate[ptm_estimate$condition=="B"&(ptm_estimate$subset=="x"),]$value <- 
                                                    fixeffects[["(Intercept)"]]
ptm_estimate[ptm_estimate$condition=="B"&(ptm_estimate$subset=="y"),]$value <- 
   fixeffects[["(Intercept)"]] + fixeffects[["subsety"]]
ptm_estimate[ptm_estimate$condition=="A"&(ptm_estimate$subset=="x"),]$value <- 
   fixeffects[["(Intercept)"]] + fixeffects[["conditionA"]]
ptm_estimate[ptm_estimate$condition=="A"&(ptm_estimate$subset=="y"),]$value <- 
   fixeffects[["(Intercept)"]] + fixeffects[["conditionA"]] + fixeffects[["conditionA:subsety"]] + fixeffects[["subsety"]]

plot_df <- rbindlist(list(pepform_df, ptm_df, ptm_estimate), fill = TRUE)
plot_points <- rbindlist(list(pepform_df, ptm_df), fill = TRUE)

plot_df$primary <- forcats::fct_inorder(plot_df$primary)

plots[[i]] <- plot_df %>% ggplot() +
  geom_line(aes(x = primary, y = value , group = rowname, color = FeatureType), size = 2) +
  geom_point(data = plot_points, aes(x = primary, y = value , group = rowname, 
                                     color = FeatureType), size = 5) +
  geom_vline(data=data.frame(x = c(43.5, 17.5, 56.5)),
             aes(xintercept=as.numeric(x)), linetype = "dashed") +
  scale_colour_manual(values = c("Peptide" = "#C3C3C3", "PTM" = "palegreen3",
                                 "PTM_estimated" = "slateblue3")) +
  labs(title = i, x = "Sample", y = "Intensity") +
  theme_light() +
  theme(axis.text.x = element_text(angle = 60, hjust=1, size = 10),
        axis.text.y = element_text(size = 16),
        legend.text=element_text(size=19),
        axis.title.y = element_text(size = 16),
        axis.title.x = element_text(size = 16),
        title = element_text(size = 19),
        strip.text = element_text(size = 16),
        legend.title =  element_blank(),
        legend.direction = "horizontal",
        legend.position = c(0.5, 0.1)) +
  annotate("text", x = 24, y = 6.5, label = "B", size = 8) +
  annotate("text", x = 70, y = 6.5, label = "A", size = 8) +
  annotate("text", x = 9, y = 6, label = "x", size = 7) +
  annotate("text", x = 30, y = 6, label = "y", size = 7) +
  annotate("text", x = 49.5, y = 6, label = "x", size = 7) +
  annotate("text", x = 73, y = 6, label = "y", size = 7) +
  ylim(-6, 6.5)

}

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

plots

## $`sp|O94769|ECM2_HUMAN (Phospho (STY)) 75`

## Warning: Removed 13 rows containing missing values (geom_point).

## 
## $`sp|O94769|ECM2_HUMAN (Oxidation (M)) 76`

## Warning: Removed 22 rows containing missing values (geom_point).

## 
## $`sp|P00747|PLMN_HUMAN (Phospho (STY)) 358`

## Warning: Removed 4 row(s) containing missing values (geom_path).

## Warning: Removed 64 rows containing missing values (geom_point).

## 
## $`sp|P01034|CYTC_HUMAN (Oxidation (M)) 40`

## Warning: Removed 29 row(s) containing missing values (geom_path).

## Warning: Removed 159 rows containing missing values (geom_point).

## 
## $`sp|P01034|CYTC_HUMAN (Phospho (STY)) 43`

## Warning: Removed 5 row(s) containing missing values (geom_path).

## Warning: Removed 40 rows containing missing values (geom_point).

## 
## $`sp|P01042|KNG1_HUMAN (Phospho (STY)) 275`

## Warning: Removed 6 row(s) containing missing values (geom_path).

## Warning: Removed 102 rows containing missing values (geom_point).

## 
## $`sp|P02765|FETUA_HUMAN (Oxidation (M)) 321`

## Warning: Removed 22 rows containing missing values (geom_point).

## 
## $`sp|P02774|VTDB_HUMAN (Phospho (STY)) 95`

## Warning: Removed 54 rows containing missing values (geom_point).

## 
## $`sp|P04075|ALDOA_HUMAN (Phospho (STY)) 39`

## Warning: Removed 1 row(s) containing missing values (geom_path).

## Warning: Removed 30 rows containing missing values (geom_point).

## 
## $`sp|P04075|ALDOA_HUMAN (Phospho (STY)) 36`

## Warning: Removed 22 rows containing missing values (geom_point).

## 
## $`sp|P05060|SCG1_HUMAN (Phospho (STY)) 311`

## Warning: Removed 1 row(s) containing missing values (geom_path).

## Warning: Removed 61 rows containing missing values (geom_point).

## 
## $`sp|P05060|SCG1_HUMAN (Phospho (STY)) 317`

## Warning: Removed 7 row(s) containing missing values (geom_path).

## Warning: Removed 97 rows containing missing values (geom_point).

## 
## $`sp|P09972|ALDOC_HUMAN (Phospho (STY)) 39`

## Warning: Removed 2 row(s) containing missing values (geom_path).

## Warning: Removed 34 rows containing missing values (geom_point).

## 
## $`sp|P10451|OSTP_HUMAN (Oxidation (M)) 284`

## Warning: Removed 37 row(s) containing missing values (geom_path).

## Warning: Removed 292 rows containing missing values (geom_point).

## 
## $`sp|P10451|OSTP_HUMAN (Phospho (STY)) 280`

## Warning: Removed 28 row(s) containing missing values (geom_path).

## Warning: Removed 269 rows containing missing values (geom_point).

## 
## $`sp|P10451|OSTP_HUMAN (Phospho (STY)) 24`

## Warning: Removed 52 row(s) containing missing values (geom_path).

## Warning: Removed 206 rows containing missing values (geom_point).

## 
## $`sp|P10451|OSTP_HUMAN (Phospho (STY)) 27`

## Warning: Removed 54 row(s) containing missing values (geom_path).

## Warning: Removed 531 rows containing missing values (geom_point).

## 
## $`sp|P10451|OSTP_HUMAN (Phospho (STY)) 195`

## Warning: Removed 102 row(s) containing missing values (geom_path).

## Warning: Removed 485 rows containing missing values (geom_point).

## 
## $`sp|P10645|CMGA_HUMAN (Phospho (STY)) 218`

## Warning: Removed 2 rows containing missing values (geom_point).

## 
## $`sp|P10645|CMGA_HUMAN (Phospho (STY)) 370`

## Warning: Removed 12 row(s) containing missing values (geom_path).

## Warning: Removed 110 rows containing missing values (geom_point).

## 
## $`sp|P10645|CMGA_HUMAN (Phospho (STY)) 161`

## Warning: Removed 17 row(s) containing missing values (geom_path).

## Warning: Removed 134 rows containing missing values (geom_point).

## 
## $`sp|P10645|CMGA_HUMAN (Phospho (STY)) 142`

## Warning: Removed 12 row(s) containing missing values (geom_path).

## Warning: Removed 72 rows containing missing values (geom_point).

## 
## $`sp|P10909|CLUS_HUMAN (Phospho (STY)) 210`

## Warning: Removed 106 rows containing missing values (geom_point).

## 
## $`sp|P13521|SCG2_HUMAN (Phospho (STY)) 106`

## Warning: Removed 23 row(s) containing missing values (geom_path).

## Warning: Removed 168 rows containing missing values (geom_point).

## 
## $`sp|P19823|ITIH2_HUMAN (Oxidation (M)) 64`

## Warning: Removed 141 row(s) containing missing values (geom_path).

## Warning: Removed 317 rows containing missing values (geom_point).

## 
## $`sp|P24592|IBP6_HUMAN (Phospho (STY)) 152`

## Warning: Removed 13 row(s) containing missing values (geom_path).

## Warning: Removed 112 rows containing missing values (geom_point).

## 
## $`sp|P24593|IBP5_HUMAN (Phospho (STY)) 124`

## Warning: Removed 30 row(s) containing missing values (geom_path).

## Warning: Removed 197 rows containing missing values (geom_point).

## 
## $`sp|P30086|PEBP1_HUMAN (Phospho (STY)) 52`

## Warning: Removed 20 rows containing missing values (geom_point).

## 
## $`sp|P51693|APLP1_HUMAN (Phospho (STY)) 515`

## Warning: Removed 14 rows containing missing values (geom_point).

## 
## $`sp|P61769|B2MG_HUMAN (Phospho (STY)) 108`

## Warning: Removed 6 rows containing missing values (geom_point).

## 
## $`sp|Q14515|SPRL1_HUMAN (Oxidation (M)) 276`

## Warning: Removed 26 row(s) containing missing values (geom_path).

## Warning: Removed 194 rows containing missing values (geom_point).

5.0.1.2 Peptidoform level

pepf_list1 <- c()
for (contrast in contrasts){
sigTable <- rowData(pe[["pepformRel"]])[[contrast]] %>% filter(adjPval<0.05)
pepf_list1 <- c(pepf_list1, rownames(sigTable))
}
pepf_list1 <- unique(pepf_list1)

plots <- list()
for (i in pepf_list1){
prot_ = rowData(pe[["pepformRel"]])[i,][["Leading.razor.protein"]]
site_ = rowData(pe[["pepformRel"]])[i,][["mod"]]
sites_ = strsplit(site_, ", ")[[1]]
for (site_ in sites_){
ptm_ = paste(prot_, site_, collapse = " ")
print(ptm_)

pepform_df <- longFormat(pe[,,"pepformRel"], rowvars = c("Leading.razor.protein", "ptm"), colvars = c("condition", "subset")) %>% as.data.frame()
pepform_df <- pepform_df %>% filter(Leading.razor.protein == prot_)
pepform_df <- pepform_df %>% filter(grepl(site_, ptm, fixed = T))
pepform_df$FeatureType <- "Peptide"
pepform_df[pepform_df$rowname==i,]$FeatureType <- "SignificantPeptide"
pepform_df <- pepform_df %>% arrange(condition, subset, rowname)

ptm_df <- longFormat(pe[,,"ptmRel"], rowvars = c("protein", "ptm"), colvars = c("condition", "subset")) %>% as.data.frame()
ptm_df <- ptm_df %>% filter(ptm == ptm_)
ptm_df$FeatureType <- "PTM"
ptm_df <- ptm_df %>% arrange(condition, subset, rowname)

ptm_estimate <- ptm_df %>% select(c("primary", "colname", "condition", "subset")) %>%
                           mutate(rowname = paste(ptm_df$rowname, "estimate"),
                                  ptm = paste(ptm_df$ptm, "estimate"),
                                  Protein = paste(ptm_df$Protein, "estimate"),
                                  FeatureType = "PTM_estimated",
                                  value = NA,
                                  assay = "model")
fixeffects <- rowData(pe[["ptmRel"]])$msqrobModels[[unique(ptm_df$ptm)]] %>% getCoef
ptm_estimate[ptm_estimate$condition=="B"&(ptm_estimate$subset=="x"),]$value <- 
                                                    fixeffects[["(Intercept)"]]
ptm_estimate[ptm_estimate$condition=="B"&(ptm_estimate$subset=="y"),]$value <- 
   fixeffects[["(Intercept)"]] + fixeffects[["subsety"]]
ptm_estimate[ptm_estimate$condition=="A"&(ptm_estimate$subset=="x"),]$value <- 
   fixeffects[["(Intercept)"]] + fixeffects[["conditionA"]]
ptm_estimate[ptm_estimate$condition=="A"&(ptm_estimate$subset=="y"),]$value <- 
   fixeffects[["(Intercept)"]] + fixeffects[["conditionA"]] + fixeffects[["conditionA:subsety"]] + fixeffects[["subsety"]]


plot_df <- rbindlist(list(pepform_df, ptm_df, ptm_estimate), fill = TRUE)
plot_points <- rbindlist(list(pepform_df, ptm_df), fill = TRUE)


plot_df$primary <- forcats::fct_inorder(plot_df$primary)

plots[[paste(i, site_)]] <- plot_df %>% ggplot() +
  geom_line(aes(x = primary, y = value , group = rowname, color = FeatureType), size =2) +
  geom_point(data = plot_points, aes(x = primary, y = value , group = rowname, color = FeatureType), size = 5) +
  geom_vline(data=data.frame(x = c(43.5, 17.5, 56.5)),
             aes(xintercept=as.numeric(x)), linetype = "dashed") +
  scale_colour_manual(values = c("Peptide" = "#C3C3C3", "PTM" = "palegreen3", 
                                 "SignificantPeptide" = "violetred1",
                                 "PTM_estimated" = "lightgoldenrod")) +
  labs(title = paste(i, ptm_), x = "Sample", y = "Intensity") +
  theme_light() +
  theme(axis.text.x = element_text(angle = 60, hjust=1, size = 10),
        axis.text.y = element_text(size = 16),
        legend.text=element_text(size=19),
        axis.title.y = element_text(size = 16),
        axis.title.x = element_text(size = 16),
        title = element_text(size = 19),
        strip.text = element_text(size = 16),
        legend.title =  element_blank(),
        legend.direction = "horizontal",
        legend.position = c(0.5, 0.1)) +
  annotate("text", x = 24, y = 6.5, label = "B", size = 8) +
  annotate("text", x = 70, y = 6.5, label = "A", size = 8) +
  annotate("text", x = 9, y = 6, label = "x", size = 7) +
  annotate("text", x = 30, y = 6, label = "y", size = 7) +
  annotate("text", x = 49.5, y = 6, label = "x", size = 7) +
  annotate("text", x = 73, y = 6, label = "y", size = 7) +
  ylim(-6, 6.5)
}
}

## [1] "sp|O94769|ECM2_HUMAN (Phospho (STY)) 75"

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## [1] "sp|O94769|ECM2_HUMAN (Oxidation (M)) 76"

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## [1] "sp|P01034|CYTC_HUMAN (Phospho (STY)) 43"

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## [1] "sp|P01034|CYTC_HUMAN (Oxidation (M)) 40"

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## [1] "sp|P01034|CYTC_HUMAN (Phospho (STY)) 43"

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## [1] "sp|P04075|ALDOA_HUMAN (Phospho (STY)) 39"

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## [1] "sp|P05060|SCG1_HUMAN (Phospho (STY)) 317"

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## [1] "sp|P10451|OSTP_HUMAN (Phospho (STY)) 280"

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## [1] "sp|P10451|OSTP_HUMAN (Phospho (STY)) 275"

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## [1] "sp|P10451|OSTP_HUMAN (Phospho (STY)) 280"

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## [1] "sp|P10451|OSTP_HUMAN (Oxidation (M)) 284"

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## [1] "sp|P10451|OSTP_HUMAN (Phospho (STY)) 280"

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## [1] "sp|P10451|OSTP_HUMAN (Oxidation (M)) 284"

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## [1] "sp|P10451|OSTP_HUMAN (Phospho (STY)) 270"

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## [1] "sp|P10451|OSTP_HUMAN (Oxidation (M)) 284"

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## [1] "sp|P10645|CMGA_HUMAN (Phospho (STY)) 142"

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## [1] "sp|P04075|ALDOA_HUMAN (Phospho (STY)) 39"

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

## Warning: 'experiments' dropped; see 'metadata'

## harmonizing input:
##   removing 360 sampleMap rows not in names(experiments)

plots

## $`_LPIVNFDYS(Phospho (STY))M(Oxidation (M))EEK_ (Phospho (STY)) 75`

## Warning: Removed 13 rows containing missing values (geom_point).

## 
## $`_LPIVNFDYS(Phospho (STY))M(Oxidation (M))EEK_ (Oxidation (M)) 76`

## Warning: Removed 22 rows containing missing values (geom_point).

## 
## $`_LVGGPMDAS(Phospho (STY))VEEEGVRR_ (Phospho (STY)) 43`

## Warning: Removed 5 row(s) containing missing values (geom_path).

## Warning: Removed 40 rows containing missing values (geom_point).

## 
## $`_LVGGPM(Oxidation (M))DAS(Phospho (STY))VEEEGVRR_ (Oxidation (M)) 40`

## Warning: Removed 29 row(s) containing missing values (geom_path).

## Warning: Removed 159 rows containing missing values (geom_point).

## 
## $`_LVGGPM(Oxidation (M))DAS(Phospho (STY))VEEEGVRR_ (Phospho (STY)) 43`

## Warning: Removed 5 row(s) containing missing values (geom_path).

## Warning: Removed 40 rows containing missing values (geom_point).

## 
## $`_GILAADESTGS(Phospho (STY))IAKR_ (Phospho (STY)) 39`

## Warning: Removed 1 row(s) containing missing values (geom_path).

## Warning: Removed 30 rows containing missing values (geom_point).

## 
## $`_GHPQEESEESNVS(Phospho (STY))MASLGEK_ (Phospho (STY)) 317`

## Warning: Removed 7 row(s) containing missing values (geom_path).

## Warning: Removed 97 rows containing missing values (geom_point).

## 
## $`_EFHSHEFHS(Phospho (STY))HEDMLVVDPK_ (Phospho (STY)) 280`

## Warning: Removed 28 row(s) containing missing values (geom_path).

## Warning: Removed 269 rows containing missing values (geom_point).

## 
## $`_EFHS(Phospho (STY))HEFHS(Phospho (STY))HEDM(Oxidation (M))LVVDPK_ (Phospho (STY)) 275`

## Warning: Removed 15 row(s) containing missing values (geom_path).

## Warning: Removed 238 rows containing missing values (geom_point).

## 
## $`_EFHS(Phospho (STY))HEFHS(Phospho (STY))HEDM(Oxidation (M))LVVDPK_ (Phospho (STY)) 280`

## Warning: Removed 28 row(s) containing missing values (geom_path).

## Warning: Removed 269 rows containing missing values (geom_point).

## 
## $`_EFHS(Phospho (STY))HEFHS(Phospho (STY))HEDM(Oxidation (M))LVVDPK_ (Oxidation (M)) 284`

## Warning: Removed 37 row(s) containing missing values (geom_path).

## Warning: Removed 292 rows containing missing values (geom_point).

## 
## $`_EFHSHEFHS(Phospho (STY))HEDM(Oxidation (M))LVVDPK_ (Phospho (STY)) 280`

## Warning: Removed 28 row(s) containing missing values (geom_path).

## Warning: Removed 269 rows containing missing values (geom_point).

## 
## $`_EFHSHEFHS(Phospho (STY))HEDM(Oxidation (M))LVVDPK_ (Oxidation (M)) 284`

## Warning: Removed 37 row(s) containing missing values (geom_path).

## Warning: Removed 292 rows containing missing values (geom_point).

## 
## $`_VS(Phospho (STY))REFHSHEFHSHEDM(Oxidation (M))LVVDPK_ (Phospho (STY)) 270`

## Warning: Removed 35 row(s) containing missing values (geom_path).

## Warning: Removed 241 rows containing missing values (geom_point).

## 
## $`_VS(Phospho (STY))REFHSHEFHSHEDM(Oxidation (M))LVVDPK_ (Oxidation (M)) 284`

## Warning: Removed 37 row(s) containing missing values (geom_path).

## Warning: Removed 292 rows containing missing values (geom_point).

## 
## $`_S(Phospho (STY))GEATDGARPQALPEPMQESK_ (Phospho (STY)) 142`

## Warning: Removed 12 row(s) containing missing values (geom_path).

## Warning: Removed 72 rows containing missing values (geom_point).

## 
## $`_GILAADESTGS(Phospho (STY))IAK_ (Phospho (STY)) 39`

## Warning: Removed 1 row(s) containing missing values (geom_path).

## Warning: Removed 30 rows containing missing values (geom_point).